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Abstract 

This paper derives a complete analytical solution for the probability distribution of the 
configuration of a non-holonomic vehicle that moves in two spatial dimensions by satisfying 
the unicycle kinematic constraints and in presence of Brownian noises. In contrast to pre¬ 
vious solutions, the one here derived holds even in the case of arbitrary linear and angular 
speed. This solution is obtained by deriving the analytical expression of any-order moment 
of the probability distribution. To the best of our knowledge, an analytical expression for 
any-order moment that holds even in the case of arbitrary linear and angular speed, has 
never been derived before. To compute these moments, a direct integration of the Langevin 
equation is carried out and each moment is expressed as a multiple integral of the determin¬ 
istic motion (i.e., the known motion that would result in absence of noise). For the special 
case when the ratio between the linear and angular speed is constant, the multiple integrals 
can be easily solved and expressed as the real or the imaginary part of suitable analytic 
functions. As an application of the derived analytical results, the paper investigates the 
diffusivity of the considered Brownian motion for constant and for arbitrary time-dependent 
linear and angular speed. 

Keywords: Brownian motion, Stochastic dynamics, Fokker-Plank equation without detailed 
balance, Stochastic processes, Localization. 
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1 Introduction 


In recent years a great interest has been devoted to investigate the statistical properties of the 
motion of active (or self-propelled) particles. These particles differ from passive particles since 
they move under the action of an internal force. Typical examples of these particles are Brownian 
motors imunnsa, biological and artificial microswimmers 0 mm mm 122 m , macroscopic 
animals min nn and even pedestrians [25). 

In the last decade, this investigation has also been addressed in the framework of mobile 
robotics, in particular by analyzing the 2D cases of the unicyle, cart, and car EDI EH, and the 
3 D case of flexible needle steering mm- In these works, analytical solutions are obtained by 
directly solving the corresponding Fokker-Planck equations by using the Fourier transform. An¬ 
other method proposed for obtaining a closed-form solution is the use of exponential coordinates 
mi eh. In all these works, the benefit of having closed-form solutions is clearly illustrated by 
introducing motion planning methods based on them. 

In this paper we want to deal with wheeled robots, which represent a special case of active 
particles moving in a two dimensional space. Their configuration is characterized by a 3D-vector 
(2 components characterize their position and one their orientation). On the other hand, in most 
of cases, a wheeled robot cannot move along any direction since it must satisfy the so-called 
non-holonomic constraints m- Our stochastic motion model (presented in section [2| differs 
from the one considered in [ l4[ EDEUE] since we assume independent errors acting on the 
linear and angular components. The model adopted in m HD [301E] refers to the case of a 
differential drive system where the errors acting on each wheel are independent. 

Starting from our stochastic motion model, we compute the analytical expression of any- 
order moment of the probability distribution (section [3]) . Specifically, each moment is expressed 
as a multiple integral of the deterministic motion performed by the mobile robot (i.e., the known 
motion that would result in absence of noise). In other words, the expression is a multiple integral 
on two time-dependent functions, which describe the time behaviour of the deterministic linear 
and angular speed. This allows obtaining any order moment for any deterministic motion, i.e., 
when the mobile robot is propelled by arbitrary time-dependent linear and angular speed. For 
the special case when the ratio between these two speeds is constant, the multiple integrals can be 
expressed as the real or the imaginary part of suitable analytic functions. In section[4]we show the 
power of the derived analytical results by investigating the diffusivity of the considered Brownian 
motion for constant and for arbitrary time-dependent linear and angular speed. Conclusions are 
provided in section [5j 

The interest in deriving the statistical properties of the motion of a wheeled robot comes 
from the possibility of improving its localization (both in precision and speed). The localization 
is a fundamental problem in mobile robotics which must be solved to autonomously and safely 
navigate. It is a non linear estimation problem. In most of cases, the localization is carried out 
by using nonlinear filters (e.g., Extended and Unscented Kalman Filters, Particle filters, etc). On 
the other hand, non linear filters are strictly connected with the solution of partial differential 
equations (see 0 [5, B] [7 and the next section [2j. Very easily speaking, for a dynamic stochastic 
process described by stochastic differential equations, it is possible to introduce a probability 
distribution which provides the probability that the process takes a given set of values. This 
probability distribution must satisfy a partial differential equation. In this sense, a non linear 
filter could be considered as a numerical solution of this partial differential equation. In the 
specific case of localization, the dynamic process is the motion of the robot together with the 
noisy data delivered by its sensors. 

We already computed the statistics up to the second order for a similar motion model HU 
US] ( in the expressions in m there are some typos, corrected in m)- Here, we extend the 
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computation in m by including any-order statistics. We also show that the expressions hold for 
a much more general model (the one described in section [2]) and not only for a simple specific 
odometry error model (as shown in [15]). As a result, this paper provides a complete analytical 
solution to a Fokker-Planck equation for which the detailed balance condition is not verifiec0 
Note that the goal of this paper is to analytically compute the moments of the probability 
distribution up to any order for the configuration of a non-holonomic vehicle. As it will be seen, 
this is a very hard task from a computational point of view. These analytical results could be 
very useful to improve the localization which is, as previously mentioned, a numerical solution 
of a partial differential equation. 


2 Stochastic motion model 

We consider a mobile robot which moves in a 2D—environment. The configuration of the mobile 
robot is characterized by its position and orientation. In Cartesian coordinates, we denote them 
by r = [x, y] T and 9, respectively. We assume that the motion of the mobile robot satisfies the 
unicycle model El. which is the most general nonholonomic model for 2D—motions, where the 
shift of the mobile robot can only occur along one direction ([cos0, sin0] T ): 


X = 

v cos 6 

y = 

usin# 

o = 

UJ 


The quantities v = v(t) and u> = ui(t) are the linear and angular speed, respectively. We assume 
that these functions are known and we call the motion that results from them the deterministic 
motion. Now we want to introduce a stochastic model that generalizes equation 0 by accounting 
a Gaussian white noise. We assume that also the noise satisfies the same nonholonomic constraint. 
The stochastic differential equations are: 

dx(t) = cos 9(t) [v(t)dt + a(t)dw r (t)\ 

dy[t) = sin 9(t) \v{t)dt + a(t)dw r (t)\ (2) 

dd(t) = w{t)dt + (3(t)dw e (t) 

where [dw r (t), dw e (t)] T is a standard Wiener process of dimension two [18]. The functions 
a(t) and /3(f) are modeled according to the following physical requirement (diffusion in the 
overdamped regime). We require that, when the mobile robot moves during the infinitesimal 
time interval dt, the shift is a random Gaussian variable, whose variance increases linearly with 
the traveled distance. This is obtained by setting a(t) = ^K r \v(t)\, with K r a positive parameter 
which characterizes our system (interaction robot-environment). Similarly, we require that the 
rotation accomplished during the same time interval dt is a random Gaussian variable, whose 
variance increases linearly with the traveled distance. Hence, we set /3(t) = y/Kg\v(t)\, with Kg 
another positive parameter which characterizes our system (interaction robot-environment). The 
equations in ([2| are the Langevin equations in the overdamped limit. From now on, we assume 
that the mobile robot can only move aheac0 i.e., v(t) > 0. We remind the reader that v(t) is a 
deterministic and known function of time. The curve length ds = v{t)dt is the curve length of the 
deterministic motion, which is independent of the Wiener process [ dw r (t ), dw e (t)] T . Hence, we 

1 Note that, deriving any-order moment, corresponds to analytically derive the probability density distribution 

m- 

2 Note that this constraint does not limit the robot motion: for instance, a go and back motion is easily obtained 
by accomplishing 180 deg rotations 
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are allowed to use the deterministic curve length ds = v(t)dt instead of the time t. By denoting 
the ratio between the angular and the linear speed by y, equations in ([2]) read: 


dx(s) = cos 9(s) [ds + Kr dw r (s)\ 
dy(s) = sin#(s) [ds + Krdw r (s)\ 
d9(s) = fj.(s)ds + Kg dw e (s) 


The associated Smoluchowski equation is: 


where: 


dp 


- V • (V lP ) + V • (X> 2 Vp) = 0 


(3) 


(4) 


• p = p(x, y, 9 ; s(t)) is the probability density for the mobile robot at the configuration 
(x, y, 9) and at time U 


• X>i = [cos0, sin0, /i] T is the drift vector; 


• T> 2 


I\ r cos 2 9 K r cos 9 sin 9 0 

4 K r cos 9 sin 9 K r sin 2 9 0 

0 0 Kg 


is the diffusion tensor. 


Our goal is to obtain the probability density p(x, y, 9 ; s). This will be done in the next two 
sections. 


3 Computation of the probability distribution 

The probability density p(x, y, 9; s) satisfies the Smoluchowski partial differential equation in 
Q, which is a special case of the Fokker-Planck equation. Since the detailed balance condition 
is not satisfied f24j, we follow a different procedure to compute p(x, y, 0; s). Specifically, we 
use the Langevin equation in © to compute the moments, up to any order, of the probability 
distribution. First of all, we remark that the third equation in © is independent of the first two, 
is independent of 9 and is linear in dw e . As a result, the probability density only in terms of 
9 , i.e., the probability Pg{9: s) = f dx J dy p(x, y, 9 ; s), is a Gaussian distribution with mean 
value 9q + f 0 ds'p(s') and variance Kgs, where 9q is the initial orientation. This same result 
could also be obtained by integrating Q in x and y and by using the divergence theorem. In 
other words, we have the following distribution for the orientation at a given value of s: 


9(s) = Af(9(s), Kgs) (5) 

where A/"(•, •) denotes the normal distribution with mean value and variance the first and the 
second argument; 9(s) = 9q + Jq ds' y{s'). Note that, according to the unicycle model, a 
trajectory is completely characterized by its starting point and by the orientation vs the curve 
length, i.e., by the function 6(s). In the following, we will call the deterministic function 9{s), 
the deterministic trajectory. 
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3.1 First and second-order statistics 

Let us consider the first equation in |3|. We compute the expression of x(s) by a direct integra¬ 
tion. We divide the interval (0, s) in N equal segments, 5s = We have: 

N 

x (s) = + e i) cosd 3 ( 6 ) 

where ej is a random Gaussian variable satisfying (ej) = 0, (ej ek) = 5jkK r 5s for j, k = 1, ■ • • , N 
(Sjk is the Kronecker delta) and 9j = 9(j5s). On the other hand, we have: 


6j = 9(j5s) + 59 m = 9j + A9j 


(7) 


m—1 


where, according to our stochastic model in ([3]), 59j is a random Gaussian variable satisfy¬ 
ing (59j) = 0, (59j ej,) = 0 and (59j 59k) = 5jkKg5s for j, k = 1,- •• ,N. As a result, 

A 9j = J2m=l m i s a ^ so a ran dom Gaussian variable satisfying (A 9j) = 0, (A 9j A 9k) = 

K( . n q £ 7j 7 ^ k q j 5 s 

K -'j j y ■ Starting from (61 and (71 and by remarking that (cos A 9j) = e ~ 

and (sin A9j) = 0, it is easy to obtain the mean value of x(s). We have: 


N 


(x(s)) = lim ^(<5s + (ej)) (cos(9,) = 

N—>oo ' 


( 8 ) 


N 


= lim > Ss cos 0j e 

N—t oo * ^ 1 


3 =1 

_ kgjSs 


f S _ Kqs' 

= / ds cos 9(s ) e 


3 = 1 


Similarly, it is possible to obtain the mean value of y(s ): 




( y(s )) = / ds' sin 9(s') e 2 


(9) 


Finally, in a similar manner, but with some more computation, we obtain all the second-order 
moments: 


( x ( s ) 2 ) 



{(1 + Xc (s')) 


cos[9(s' + s") — 0(s 7 )] — Xs( s ') sin[0(s' + s") — 0(s 7 )]} + 



{y(s) 2 )=j ds' [ ds”e {(1 — Xc( s ')) 

Jo Jo 

cos[9(s' + s”) - 0(s 7 )] + Xs(s') sin[0(s 7 + s") - 0(s 7 )]} + 



(*(«) y(s)) 



{x*W) 


( 10 ) 


(ii) 


( 12 ) 
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cos[0(s' + a") - 0(s')] + Xc(s') sin[0(s' + s") - 0(s')]} + 
K r 


+ - 


ds' Xs(s') 


<7x0(s) = (x(s) 8(s)) - (x(s)) 8(s) = 2K e 


<T y g{s) = ( y(s) 8(s)) - ( y{s )) 8(s ) = -2 K e 


d(y(s)) 

0K e 

d(x(s)) 


dKg 


(13) 

(14) 


where Xc{s') = cos[20(s 7 )] e~ 2Kes and Xs(s') = sin[20(s 7 )] e~ 2Kes . 

Obtaining the expression of higher-order moments will demand more tricky computation and 
will be dealt separately, in the next subsection. Here we conclude by considering the quantity 
D(s) 2 = x(s) 2 + y(s) 2 , which provides the time-evolution of the square of the distance of the 
mobile robot from its initial position. From (10) and (111 we obtain a simple expression for its 
mean value: 


K r s + 2 ds' 


(D(s) 2 ) = 

ds costs' + s") - 0(s')] 


(15) 


3.2 Computation of any-order moment 

We introduce the following two complex random quantities: 

N N 


0) = JV 1 ™.E^ s + e i) e ^ ; W ( S ) = E 1 ^ + e l) e 


- i6n 


3 =1 


N—>oo 1 


(16) 


3 = 1 


From Wo\ and (161 it is immediate to realize that a;(s) = u ( s )+ w ( s ) _ Similarly we have y(s) = 
u(s)- w[s) ' jj ence> j n orc (er to compute any-order moment which involves x(s) and y(s) it suffices 
to compute (u(s) p w(s) q ) for any p,q £ A f. The computation of this quantity requires several 
tricky steps, which are provided in appendixjA] The key is to separate all the independent random 
quantities in order to compute their mean values. This is obtained by arranging all the sums in 
a suitable manner (see appendix [A~| for all the details). We have: 


2 n 


L4 q -i 

l ( s yw( s y)= 


p 


2n — l J \ l 


(17) 


E 


2 n — l\ { l 


m 


m 


n=0 1=0 

m!(2n — l — to — 1)!!(Z — m — 1)!! s ri 


i(p-q)8 0 


2n — l — m 


l — TO 


\(p — 2n + l)\(q — l)\ 


E 


ds i / ds 2 
J Si 


rs r/ 3 - 1 

• / dspexp E [i(p - q + $ b )- 
dsp-i lb=0 


■[0(s 6+ i) - 0(s 6 )] - 


(p-q + < Fb) 2 (s b+ 1 - S b )Kg 


where: 
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• L^J is the largest integer not greater than 7 s ; 

• the second sum on l (i.e., Ya=o) i s restricted to the values of l for which 2 n — l < p and 
l < q (or, equivalently, we are using the convention that (®) = 0 when y > x)\ 

• the sum on m (i.e., Yh-m) g° es from 0 to the minimum between l and 2 n — l and it is 
restricted to the integers m with the same parity of l (hence, both 2n y(~ m and are 
integers); 

• the symbols "!" and "!!" denote the factorial and the double factorial, respectively [5] (note 
that 0! = OH = (-1)!! = 1); 

• (3 = p + q — n — mis the dimension of the remaining multiple integral (note that a multiple 
integral of dimension m has already been computed and provided the term s m ); 

• the sum on c, i.e., Y2& * s the sum over all the vectors c of dimension j3, whose entries 
are —2, —1, 1 and 2: specifically, each vector c has l -^ rk entries equal to 2, q — l equal 
to 1, p — 2n +1 equal to —1 and 2ra ~^~ m equal to —2 (note that the sum Ylc consists of 
(/,) (/-?«,) (1?) “Mends); 

• <f>6 = J2a=l ( note that <&p = q-p)\ 

• Sq = 0 and 9(sq) = 0q (note that sq is not a variable of integration). 


In order to complete the derivation of the statistics for our problem, we need to compute any- 
order moment which also involves the orientation 9. We provide the formula for the quantity 
(u(s) p w(s) q 9(s ) r where 9(s) = 9(s) — 9(s). The details of this computation are provided in 
appendix |B| We have, for any p,q,r £ Af\ 


L^J 2r, 


l{s) P w(s) q 9{s) r ) = Y. K rY 


n—0 1=0 


2 n — l) \ l 


Y (^ n m Qjm!(2n -l-m- !)!!(/ - m - 1)!! s m e ^ p - q)0 ° 


2n — l — m \ ! (l — to 


!(p-2n + 0!(9-0!EE 


r!(7/3+i - 


c 7 

^73+1 


nLo 7b+ 


11/(2 ps ps PS 

—— / dsi / ds '2 ■ ■ ■ / dsp(s — sp) — 
l! J 0 J si J sp-\ 


exp < iY(P ~ 1 + $ b )[9(s b+ 1 ) -0(s 6 )] 


b= 0 


3-1 


6=0 
| '?b+ 



(p-q + $ b ) 2 Ke(s b +i-s b y 

S *6/ t-Xp 

2 


7 b+i! (i{p- q + ^b)V K e(s b + 1 - s&)) 

' -J 


7fe+i —2a ' 


a—0 


a!(7 6+ i - 2a)!2“ 


(18) 












where the sum over 7 , i.e. ^71 the sum over all the vectors 7 = [71,72, • ■ ■ , 7^, 7/3+1], where 7 g 
(b = 1, • • ■ , ft) are positive integers and 7 / 3+1 is a positive integer with even parity. Additionally, 
they satisfy the constraint: X^ 6 =i 7 & = r - 


4 Diffusivity for constant and for arbitrary time-dependent 
linear and angular speed 

In this section we want to illustrate the power of the formulas derived in the previous section, 
which hold for any time-dependent linear and angular speed, to compute the statistics in several 
specific cases. In particular, we compare the results obtained by using our formulas with the 
results that it is possible to obtain by running Monte Carlo simulations. We generate many trials 
of random motion via Monte Carlo simulations according to the motion model in <©• The initial 
configuration is the same for all the trials and it is a;(0) = y(0) = 0(0) = 0. The final curve 
length is also the same and it is s = 1. We consider several values for the two parameters Kg and 


K r and we set /z(s) constant in 4.1 and variable in 4.2 We will see that, as the number of trials 


increases, the statistical properties of the motion directly obtained from the trials approach the 
ones obtained by using our formulas. In particular, we will refer to the motion diffusivity and we 
compute the moments up to the forth order. We start by computing for each trial, the square of 
the final distance from the origin, i.e. the distance at s = 1 : 


D{l) 2 = x{lf + y{lf 


(19) 


Once we have the previous quantity for many trials, we compute its mean value and its variance. 
On the other hand, our formulas allow us to directly obtain both the mean value and the 
variance. Specifically, equation (15 ) provides the mean value. Regarding the variance we have: 
ctd ( s )2 = (D(s) 4 ) — (-D(s) 2 )“. Additionally, D(s ) 4 = :r(s ) 4 + 2 x(s) 2 y(s) 2 + y(s) 4 = u(s) 2 w(s) 2 , 
where we used cc(s) = U ( B )+ W ( S ) an d y( s ) = u ( s )~. w ( s ) . qq ie q uan tity ( u(s) 2 w(s ) 2 ) can be easily 
computed by using equation © with p = q = 2 (see appendix [C| for the details). Its analytical 
expression includes the following six terms: n = 0 , Z = 0 , ?n = 0 (C*ooo); n=l,Z = 0 ,m = 0 
(C 100 ); n = 1, l = 1, m = 1 (Cm); n = 1, l = 2, m = 0 (C 120 ); n = 2, l = 2, m = 0 (C 220 ) and 
n = 2, l = 2, m = 2 ( 0222 )- We have: 


(D(s) 4 ) — 0*000 + 0*100 + C *111 + C*i20 + 0*220 + 0*222 ( 20 ) 

See appendix [C] for the analytical expression of the previous terms. In the next two subsections, 
we compare the results for the mean value and the variance of D( l ) 2 obtained by running many 
trials and by using our analytic expressions. In |4.1| we refer to a constant ratio between the linear 
and the angular speed while in |4.2| we refer to a generic varying ratio. Note that this analysis is 
a test for the 4 th order statistics of the considered Brownian motion. 


4.1 Constant ratio between linear and angular speed 

When /Lt(s) = /tq we have that 0(s) is linear in s: 


9(s) = 9 0 + y 0 s 


( 21 ) 


and, the corresponding active trajectory, is a circumference with radius ^ (note that when the 
angular speed vanishes, the radius becomes infinite and the active trajectory becomes a straight 
line). The computation of the integrals in (|T7|) and ( 181 is immediate. As a result, the expression 
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of <D( S ) 2 > can be easily provided in terms of fj,o, s, K, and Kg. Let us introduce the following 
complex quantity: 


- Ks ^- 
z = +*/Li 0 


By a direct analytical computation of the double integral in (151 we easily obtain: 


( 22 ) 


( D(s ) 2 ) = K r s + 2SR j j, | (23) 

where the symbol 3?{-} denotes the real part of a given complex quantity. 

We do not provide here the analytical expression for ^Z?(s) 4 ). We only remark that all the 
integrals appearing in the expressions in appendix [C] can be computed in a similar manner and 
are expressed as the real part of suitable analytic functions of the following complex quantities: 
z, + ifJ-o and —2 Kg + 2i/Li 0 - 


E 



0.4 

0.3 

0.2 

0.1 

0 


-0.3 -0. 


2 - 0 . 


Figure 1: Trajectories generated in the case of constant ratio between linear and angular speed. 
The blue line is the trajectory generated without noise (i.e, Kg = K r = 0) and the other lines 
are five trajectories obtained by setting Kg = K r = 0.01 to. 


In the following, we illustrate some numerical results obtained by setting /r(s) = H o = 5. 
We consider two cases of Brownian noises, the former is characterized by Kg = K r = O.Olm 
and the latter by Kg = K r = 1 m. Figure |T] displays the trajectory generated without noise 
(blue line) together with five trajectories obtained by setting Kg = K r = O.Olm. For the 
trajectory without noise, the final -D(l) 2 is equal to 0.0573 m 2 . We use the expression in (231 
to compute the mean value and a similar expression, which depends on the three mentioned 
complex quantities z, + i/j.o and —2 Kg + 2i/iQ, to compute the variance. We obtain the 
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II 

, ® 

trials 

{D(lf) [m z ) 

a D(i ) 2 

0.01 

10 3 

0.0688 

0.0011 

0.01 

10 4 

0.0664 

0.0012 

0.01 

10 5 

0.0679 

0.0012 

1 

10 3 

1.1189 

1.3736 

1 

10 4 

1.1124 

1.2767 

1 

10 5 

1.1126 

1.3035 


Table 1: Constant ratio between linear and angular speed. Values of (D( l) 2 ) and cr 2 D (yy obtained 
by running 10 3 , 10 4 and 10 5 trials. 


following values: (Z3(l) 2 ) = 0.0680 m 2 and cr^ny = 0.0012 m 4 when Kg = K r = 0.01m and 
(D(l) 2 ) = 1.1130 m 2 and = 1.3052 m 4 when Kg = K r = lm. 

We ran 10 5 trials and we computed from them the mean value and the variance. Figure [ 2 ] 
displays the values of D{ l) 2 obtained in the first 10 4 trials together with the value of -D(l) 2 for 
the trajectory without noise (red line), the value of (Z3(l) 2 ) obtained by averaging D( l) 2 on 
the 10 4 trials (green line) and the value of (U(l) 2 ) obtained by our formulas (black line). The 
image on the left refers to the case Kg = I\ r = 0.01m and the image on the right to the case 
I\g = K r = lm. 



trial trial 

Figure 2: Constant ratio between linear and angular speed. Values of D( l) 2 for 10 4 trials together 
with the value of Z?(l) 2 for the trajectory without noise (red line), the value of (D( l) 2 ) obtained 
by averaging D( l) 2 on the trials (green line) and the value of (Z?(l) 2 ) obtained by our formulas 
(black line). Kg = K r = 0.01m and Kg = K r = lm in the left and right image, respectively. 


Table |Tj reports the values of (H(l) 2 ) and (J 2 D yy obtained by running 10 3 , 10 4 and 10 5 trials. 
It is possible to see that, as the number of trials increases, the results converge to the values 
obtained by using our formulas. 
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Figure 3: Trajectories generated in the case of varying ratio between linear and angular speed. 
The blue line is the trajectory generated without noise (i.e, Kg = K r = 0) and the other lines 
are five trajectories obtained by setting Kg = K r = 0.01?n. 
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4.2 Arbitrary time-dependent ratio between linear and angular speed 

We consider now a generic function /i(s). Specifically, we consider several kinds of dependence 
on s obtaining very similar results. In the following, we illustrate the results obtained by setting 
/r(s) = 10 s. As in the previous case, we consider the two cases of Brownian noises, i.e., charac¬ 
terized by Kg = K r = 0.01m and by Kg = K r = lm. Figure [3] displays the trajectory generated 
without noise (blue line) together with five trajectories obtained by setting Kg = K r = 0.01m. 
For the trajectory without noise, the final D( l ) 2 is equal to 0.1021 m 2 . In this case we directly 
use the expressions in (151 in (20 1 and in appendix [C| by numerically computing all the inte¬ 
grals. We obtain the following values: (D(l) 2 ) = 0.1124 m 2 and = 0.0026 m 4 when 

Kg = K r = 0.01 m and (D(l) 2 ) = 1.1443 m 2 and cr^^y = 1.4681 m 4 when Kg = K r = lm. 

We ran 10 5 trials and we computed from them the mean value and the variance. Figure [4] 
displays the values of -D(l) 2 obtained in the first 10 4 trials together with the value of -D(l) 2 for 
the trajectory without noise (red line), the value of (_D(l) 2 ) obtained by averaging D( l) 2 on 
the 10 4 trials (green line) and the value of (D(l) 2 ) obtained by our formulas (black line). The 
image on the left refers to the case Kg = K r = 0.01m and the image on the right to the case 
Kg = I\ r = lm. 



trial trial 


Figure 4: Varying ratio between linear and angular speed. Values of D( l ) 2 for 10 4 trials together 
with the value of Z?(l ) 2 for the trajectory without noise (red line), the value of (l?(l) 2 ) obtained 
by averaging D{1) 2 on the trials (green line) and the value of (.D(l) 2 ) obtained by our formulas 
(black line). Kg = K r = 0.01m and Kg = K r = lm in the left and right image, respectively. 


Table |2j reports the values of (Z?(l) 2 ) and obtained by running 10 3 , 10 4 and 10 5 trials. As 

in the case of constant ratio between linear and angular speed, we remark that, as the number 
of trials increases, the results converge to the values obtained by using our formulas. 
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II 

, ® 
S*5 

trials 

(D(l)^) (m z ) 

a D(iy (™ 4 ) 

0.01 

10 3 

0.1139 

0.0022 

0.01 

10 4 

0.1107 

0.0026 

0.01 

10 5 

0.1123 

0.0026 

1 

10 3 

1.1094 

1.2534 

1 

10 4 

1.1605 

1.5289 

1 

10 5 

1.1494 

1.4708 


Table 2: Varying ratio between linear and angular speed. Values of (D(l) 2 ) and obtained 

by running 10 3 , 10 4 and 10 5 trials. 

5 Conclusion 

In this paper we derived the statistics, up to any order, for 2D—Brownian unicycle dynamics. 
The chosen kinematic constraint is modelled by the unicycle differential equation which is a very 
general constraint for a 2D motion. According to this model, the mobile robot can freely rotates 
but only one direction for the shift is allowed. This model is suitable to characterize the dynamics 
of many wheeled robots: namely, the ones that satisfy the unicycle dynamics. Additionally, the 
considered deterministic motion is very general. The expressions here provided for the statistics 
hold for any deterministic trajectory that satisfies the mentioned kinematic constraint. The 
expressions contain multiple integrals over the deterministic trajectory. 

We showed the power of the derived analytical results by investigating the diffusivity of the 
considered Brownian motion for constant and arbitrary time-dependent linear and angular speed. 

We want to remark that this paper provides the analytical expressions for the statistics, up 
to any order, of a non-trivial Brownian motion, i.e. propelled by arbitrary force and torque. To 
the best of our knowledge, analytical solutions for any order moment in the case of arbitrary 
linear and angular speed have never been derived in the past. 
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A Computation of ( u(s) p w(s) q ) 

We have: 

( u(s) p w{s) q ) = lim Y e <*q')\ (24) 

N^-oo 2 —' \ / 

jl ’“jpk 1 ‘“kq 

({6s + e jl )--- {6s + ej p ) {5s + e kl ) ■ ■ ■ {5s + e kq )) 
where each index goes from 1 to N. Let us focus our attention on the quantity 


((<5s + e h ) ■ ■ ■ {6s + e jp ){5s + e kl )---{6s + e kq )) 

We can write this quantity as the sum ofp + q + 1 terms, i.e., 5s p+q Co + 5s p+q ~ 1 C\ + 5s p+q ~ 2 C2 + 
• • • + 6sC p+q -1 + C p+q . We trivially have Co = 1. Concerning the remaining coefficients, we 
remark, first of all, that only the ones with even index are different from 0. In other words, 
C 2 n+i = 0, n = 0,1, • ■ • , L^J (where pj is the largest integer not greater than x). Let us 
compute C 2 n- This coefficient is the sum of ( p p 9 ) elements, each of them being the average 


of a product of 2n terms "e". On the other hand, since the exponential in (24) is symmetric 
with respect to the change j a f> )„■, a, a' = 1, • • ■ ,p and with respect to the change kj, -f-> hy , 
b, b' = 1, ■ • • , q but it is not symmetric with respect to the change j a -e-t kb, we need to separate 
the ( P 2 n) elements in several groups. Specifically, let us denote by l the number of e with index 
of type k. The number of elements belonging to this group is ( 2 p _ ; ) ( f). Note that we set (£) = 0 

when b > a. Note also that E?=o ( 2 „-j) (?) = ( P 2 n )■ We now remark that, because of the 
statistical properties of e, the average of the product of 2 n terms e is different from zero only 
when their indexes are equal two-by-two. Hence, we have to consider all the combinations of 
products of terms e, whose indexes are equal two-by-two and which differ for at least one pair 
of indexes. On the other hand, for a given element in the group characterized by a given n and 


l, the effect in (24) depends on the number of pairs which are hetero (i.e., with an index of type 
j and one of type k). Let us denote by m the number of pairs which are hetero. Note that m 
has the same parity of l. Indeed, by definition we have l indexes of type k. Additionally, we are 
using to indexes of type j and to indexes of type k to make to hetero pairs. Hence, l — m indexes 
of type k remain and, with them, we have to make homo pairs of type k. Similarly, we have 
to make 2ra ^~ m homo pairs of type j. Hence, must be integer. 

We compute the number of combinations of the products of 2 n terms e, with l indexes of 
type k, where the indexes are equal two-by-two, with to pairs which are hetero and which differ 
for at least one pair. We obtain: ( 2 ? p z ) {^ n )m\{2n — l — to — 1)H(Z — to — 1)!!. Following this, we 
can write (241 as follows: 


( u{s) p w{s) q ) = lim Y y 5s p+q 

N —>00 LL ' 

{j} F {k} q "=0 

2 n / \ / \ min(l,2n—l) 

e j,) :) e 


—2 n 


2n — l\ ( l 


. 2 n — lJ \lJ * —' \ to 

/—0 m —0 odd/even 

(2 n-l-m- 1 )!!(Z - to - 1 )!! ( e 2 h ) 6 jlh (e? 3 ) 6 j3j4 ■ ■ ■ 

e j 2n -l-m.-l ) ( C fcl ) ^*1*2 ( C fc 3 ) ^*3*4 


(4_m +1 ) 5 *l 


-m + lj2n-l-m + l 
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(25) 







<y<»«-»■•>) 


where, for the brevity sake, we denoted by Y^{j} {k} the sum k t ---k =v Each average (e 2 ) 

provides K r Ss. Hence, all together, they provide K™5s n . Additionally, the number of Kronecker 
deltas is n. Hence, in the limit of N —> oo for each value of n we get a multiple integral of 
dimension p + q — n. Note that, when the indexes are equal h-by-h (h > 2), the result in the 
limit N —^ oc vanishes since, the power of 5s, is larger than the number of sums. By a direct 


computation in (251 we obtain: 


(u(s) p w(s) q ) = 


L2 n 

E a '.”E 

n =0 1=0 


P 

2 n — l 


min(l,2n—l) 

E 

m—0 odd/even 


(26) 


2 n — l 
m 


m\(2n — l — m — 1 )!!(/ — m — 1 )!! s m lim 

N—>oo 


E 

{j°}A3 d }p{k°} x {k d }r 


5s p 


exp 


{*[% 


+ OjB A 


2 ( 0 ,, 


"+ + Qj*) ~ ^ — 2(0 fc , + • • • + 0 fc d)]|^> 


2 n—l—m 


= p— (2 n - l ), ?? = l -~Y k , X = Q^z) and /3 = p + a + 7j + x= P + q- n — 


where p = 

We must compute the average of the exponential in (26) and then we compute the limit N 


m. 

oo. 

We remark that the various 6 in the exponential contain random quantities (i.e., the 59 at 
different time steps). In order to proceed we have to separate all the random quantities which are 
independent. We start this separation by redefining the indexes in the sum { j ,} {k s } {k d } • 

Specifically, we consider the new indexes j S j d k k which differ from j s j d k s k d since they are 
ordered (in increasing order). For instance, j 1 < j 2 < ■ ■ ■ < j a . Hence, the last sum in (26 1 


can be replaced with P !ct! 0 ! X ! The f ° Ur typeS ° f 

indexes are not ordered among them. Hence, the sum includes all the possible combinations which 
maintain the order only restricted to a single index type. For instance, a possible combination 
is: k 1 < k 2 < j * < k\ < j 2 < k 3 < < j 2 < • • •. We introduce the /3 ordered indexes 


w\ < u >2 < • • • < wp. Additionally, let us denote with = J2c=a +1 ^c an d — E c =a+i c = 
9(b5s) — 9(a5s), where 0(s) is the deterministic trajectory. Finally, let us define a = 2p+a —x~2rj. 


The sum in the exponential in (261 contains a(A q 1 + A 0 1 ). Then, depending on the combination 
of the indexes j s j d k s k d , we have a different result for the (A(fJ +A^). Specifically, if w\ = jf the 
sum in the exponential contains (a — 2)(AA^ + A^). If vJi = it contains (a — 1)(A^ + A^). 

If Wi = fc* it contains (a + 1)(AAJ + A^). If Wi = kf it contains (a + 2)(AA^ + A^). We 
introduce the vector c = \c\, ■ ■ ■ , cp] T whose entries are —2, —1, 1 and 2. Specifically, it contains 
p entries equal to 2, x equal to 1, a equal to —1 and p equal to —2. It is easy to realize that we 
have ( /3 ) (^” x ) ( P+I? ) vectors c. Finally, we define <!>{, = X)a=i c a- Note that = — a. According 


to this, the last sum in (261 can be written as follows: 


e'^-^plc r!??!x!^ ^ 5s?- 


JaS o , 


exp |* E [(“ + <J> b)(‘^ 
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Hence, we have: 


L^J 2u 


i{sYw{sY)= Y. K rY 


n —0 /—0 


p q 
2n — l) \l 


(27) 


min(l,2n—l) 

E 

m =0 odd/even 


2n — l\f l 
m J \m 


i!(2 n — l — m — 1)!!(Z — m — !)!!• 


■s m e iae ° 


p\a\p\x^^Y E 


' N —>oo 

c 


exp | * E [(“ + $ b)( A ^ +1 + A S!) +1 ) 


b=0 


Now we can compute the average since we were able to separate all the independent quantities. 
We have: 


expliY [(« + $ b)( A ^ +1 + A ®E' 


6=0 


exp | * Y [( a + EE 

l 6=0 

II (exp{* [(a + $&)A“‘ +1 ] J 
' 0 

exp\i Y [(« + $ b ) A ^ +1 


/ 3-1 


6=0 


6=0 


n exp < - 


6=0 


(a + $ b ) 2 (iu b+ i - w b )K e 5s 


We used the equality (77 = Af (0,1)): 


(**”) =' 


(28) 


(29) 


with A = *(a: + <h b )cr b+ 1 and cr^ +1 = (w b +i — w b )Kg5s. By substituting equation (28) in (271 and 
by taking the limit (TV —> 00 ), we finally obtain: 


L4H 2n 


p \fq 
2 n — l) \l 


min(l,2n—l) 

• E 

m =0 odd/even 


(s) P w{s) q ) = Y K rY 

=0 ;=0 

m\(2n — l — m — !)!!(/ — m — !)!!• 


( 30 ) 


2 ri — Z\ / l 
m J \m 
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s"V“% 


Mr/'-X'- E / ds i [ ds 2 - - f dsp 

0 J s i J s a i 


i(a+4>6)[fl(s6 + i)-e(s6)]- 


(o + f b ) 2 ( ab+1 - ab )f 


which coincides with ((171). 


B Computation of (u{s) p w(s) q 9 r ^ 


This computation follows the same initial steps carried out in appendix [X] We 
pression equal to the one given in (271 but the mean value at the end must be 

{*Ef=o [(«+^)(A“r+ a^: + i )]} 


exp < 


9 r ). The deterministic part can be 


the mean value. We need to calculate: 


obtain an ex- 
replaced with 
factorized out 


exp < i 


E [o+ $ b) A - 


Wb+1 

Wb 


6=0 

3-1 


9 r ) = 


dr n { exp *( a+$ b) A ®r i } 


76+1 


6=0 

By using the multinomial theorem we can write: 

( 0 _ V 0 

Mixr =i>n 

\b =0 / 7 6=0 

where uip+i = N and the sum over 7 is the sum over all the vectors 7 = [71, • • • , 7^+1], where 7 b 
(&=!,••• , /3 + 1) are positive integers satisfying the constraint: ^ 6=1 76 = J". Hence, we obtain 


(a;:* 1 ) 


76+1! 


3-1 


<?n 

\ 6=0 


i(a+3> b )A™ b + 1 


(31) 


E r! 

7 


( A g + 0 

7/3+1! 

= E 


n 

6=0 

r! 


i(a+$ 6 )A 

e b 


( A -6+l\ 

-1 y w b J 


76 + 1 


76+1! 


(Asr 1 ) ) 


7 nf =0 76+i! 

( A g +i ) 7b+1 


6=0 


We use the following two standard results for a normal distribution: 

0 7/3+1 °dd 


(a m 


a p+i (7/9-J-i — !)!! 7 / 3+1 euen 


(32) 


and 


exp 


^b+i 


i(a + d> b )A- 


( A Sr) 7b+1 ) 


(33) 
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= <?b+l ex p 


(a + $h) 2 crfc +1 


E 

a =0 


76+i ! ( »(g + $b)vb+i) lb+1 2a 

+76+1 - 2a)!2 a 


with fjj , 1 = /vgJs(?z;t, + i —w;,). Equation (331 is obtained starting from (291 and by differentiating 
7 {, + i times with respect to A. Hence, the expression in (311, becomes: 


9 r jezp i(a + $ b )A^ +1 J } 


(34) 


6=0 


= E 

7 


(7/3+1 - 1)!! W f ^ _wS±i 


nLo76+i! 


nw b + ?*- : 
6=0 


L^J 

E 

a =0 


76 + i! (»(q + 4» b )g fe+1 ) 
«!(76+i - 2a)!2° 


7!)+l-2a 


where now the sum over 7 only includes the vectors 7 whose last entry is even. By taking the 
limit N —» 00 , this expression becomes: 


£ 




7 Il6=o76+i! 


3-1 


6=0 


1—r I , , Ti>+1 (“+*b) -K+Ob+l-Sb) 

11 \ (Sb+i - Sh) * e~ » ■< 


76 +i! (i{at + <^ b )\/Ke(sb+i - s b /j 


7b+i —2a - 


a —0 


and we finally obtain: 


a!(76+i ^ 2a)!2 a 


L^srM 2n 


*(«) p w(sye( s y) = £ 


n=0 Z=0 


p 


min(l,2n—l) 

E 

m =0 odd/even 


2n — l\ f l 
m ) \m 


■^“V^WEE 


2n — l) \l y 
m\(2n — l — m — 1)!!(Z — m — 1)!! 
7 "-(7/3+i - l )'-'- K e 


(35) 


ds 1 / ciso- 


IK= 0 W 7o 

dS/3 ( s _ 


c 7 1 16=0 

2/3+1 


' S/3-1 
0-1 


6=0 


1- I x T'b+i (<x+®b) K e( s b +1 ~ s b) 

ilUs6+i-S6) 2 e 
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J 7b+l! (i{Oi + $b)\/Kg(Sb +1 - s b )^j 

' -J 

a—0 


76 + i-2q ' 


which coincides with (181. 


a!( 7 h+ i - 2 a)!2 a 


C Computation of (^.D(s) 4 ) = (u(s) 2 w(s)‘ 


We use equation (171 with p = q = 2 to obtain the analytical expression of (_D(s) 4 ) = 
(u(s) 2 iu(s) 2 ). The first sum on n includes the three terms n = 0,1,2. When n = 0, the 
second sum on l and the third sum on to only include one term, i.e., l = m = 0. All the factors 
before the sum on c are 1 with the exception of (p — 2n + l)\ = 2 and (q — l)\ = 2. Additionally, 
/3 = 4 and the sum on c consists of six vectors, i.e., [—1,-1,1,1], [—1,1,—1,1], [—1,1,1,-1], 
[1, —1, —1,1], [1, —1,1, —1], [1,1, —1, —1], Hence we obtain the following term: 


Cooo — 8 


ds\ 


dso 


dsz 


ds4 |( 


fi: e(~ s l~ 3s 2+ 3s 3 + s 4) 


COS (— 6 \ — 82 + 83 + 84 ) + 


+e 2 [cos (— 9\ +02 — 03 + 04 ) + cos (—0i + 02 + 03 — 04 )] ^ 

where Qj = 6(sj), j = 1, 2, 3,4. 

Let us consider now the terms with n = 1. In this case, l = 0,1,2. For each value of l we 
have in general several values of to, which must have the same parity of l. Hence, when 1 = 1, 
we only have to = 1. When l = 2, we can have m = 0, 2. On the other hand, when l = m = 2, 
the factor ( 2? = (°) vanishes. By an explicit computation, it is possible to see that C 120 is 
the conjugate of Cioo- Hence, their sum is real and it is: 


C 


100 


+ C120 — 4A r / ds\ / ds2 / ds3 j c 
J 0 >/ Si J So. 


Kq ( —4s 1 +3s2+ s 3) 


COS (201 — 0 2 — 03 ) + 


+e 


Kg(-si+S3) 


cos ^ 2$2 — ^ 3 ) H - c 


Kq( — s 1 — 3s2+ 4s 3) 


cos (—0i — 02 + 20; 


0 } 


The last term with n = 1 is the one with l = m = 1. By a direct computation we obtain: 


C 111 = 8 K r 


ds i 


ds2e 


Kq(-s 1 +s 2 ) 


COS (-01 + 02) 


Finally, when n = 2 the sum on l includes the five values: l = 0,1, 2,3,4. On the other hand, 
only the term l = 2 does not vanish since the product ( 2 p _,) (?) is zero in all the other cases. 
When l = 2 we have two possible values for to, i.e., m = 0, 2. We obtain the following two 
contributions: 



C222 = 2 AT 2 s 2 
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